clear; clc; close all;

%% 参数设置
k = 3;
c1 = 1e-9; c2 = 0; c3 = 0;

mu_list = linspace(-3*pi, 3*pi, 50);
eta_list = linspace(-3*pi, 3*pi, 50);

h = 0.01;
T = 100;
transient = 50;
N = T / h;
N_trans = transient / h;


mu_list = linspace(-3*pi, 3*pi, 100);
eta_list = linspace(-3*pi, 3*pi, 100);

LE_max_vec = zeros(length(mu_list)*length(eta_list),1);

%% 并行扫描mu-eta空间，计算最大李雅普诺夫指数
parfor idx = 1:length(mu_list)*length(eta_list)
    [i, j] = ind2sub([length(eta_list), length(mu_list)], idx);
    mu = mu_list(j);
    eta = eta_list(i);
    
    x0 = [0; 0; 0];
    LEs = computeLEs(k, c1, c2, c3, mu, eta, h, N, N_trans, x0);
    LE_max_vec(idx) = LEs(1);
end

LE_max_map = reshape(LE_max_vec, length(eta_list), length(mu_list));

%% 绘制热图
figure;
imagesc(mu_list, eta_list, LE_max_map);
set(gca,'YDir','normal');

n_neg = 50; 

% 正值区颜色数量
n_pos = 50;

% 创建负值区黑色
neg_colors = zeros(n_neg, 3); % RGB全0

% 创建正值区渐变色（比如黄色到红色）
pos_colors = [linspace(1,1,n_pos)', linspace(1,0,n_pos)', linspace(0,0,n_pos)'];

% 合并colormap
custom_cmap = [neg_colors; pos_colors];

% 应用colormap
colormap(custom_cmap);

% 限制颜色轴范围
caxis([-0.25 0.25]);

% 显示colorbar
colorbar;

caxis([-0.25 0.25]);  % 论文色条范围
xlabel('\mu');
ylabel('\eta');
title('Fig.7b: Maximum LE in \mu-\eta plane');
xticks(-3*pi : pi : 3*pi);
xticklabels({'-3\pi','-2\pi','-\pi','0','\pi','2\pi','3\pi'});
yticks(-3*pi : pi : 3*pi);
yticklabels({'-3\pi','-2\pi','-\pi','0','\pi','2\pi','3\pi'});

%% --- 李雅普诺夫指数计算函数 ---
function LEs = computeLEs(k, c1, c2, c3, mu, eta, h, N, N_trans, x0)
    n = length(x0);
    X = x0;
    Q = eye(n);
    sum_LE = zeros(n,1);
    
    for step = 1:N
        X = RK4_step(@(x) reduced_system_ode(x,k,c1,c2,c3,mu,eta), X, h);
        J = jacobian_matrix(X,k,mu,eta);
        Q = Q + h*J*Q;
        [Q,R] = qr(Q);
        
        if step > N_trans
            sum_LE = sum_LE + log(abs(diag(R)));
        end
    end
    
    LEs = sum_LE / (h*(N - N_trans));
    LEs = sort(LEs, 'descend')';
end

%% --- 降维模型系统方程 ---
function dx = reduced_system_ode(x,k,c1,c2,c3,mu,eta)
    X = x(1); Y = x(2); Z = x(3);
    dx = zeros(3,1);
    
    dx(1) = Y + Z - k * sin(Y + eta) + k * sin(eta) + c1;
    dx(2) = -X + Z + c2;
    dx(3) = -X - Z + k * sin(X + mu) - k * sin(mu) + c3;
end

%% --- 雅可比矩阵 ---
function J = jacobian_matrix(x,k,mu,eta)
    X = x(1); Y = x(2); Z = x(3);
    J = zeros(3,3);
    
    J(1,1) = 0;
    J(1,2) = 1 - k * cos(Y + eta);
    J(1,3) = 1;
    
    J(2,1) = -1;
    J(2,2) = 0;
    J(2,3) = 1;
    
    J(3,1) = -1 + k * cos(X + mu);
    J(3,2) = 0;
    J(3,3) = -1;
end

%% --- RK4一步积分 ---
function x_next = RK4_step(f,x,h)
    k1 = f(x);
    k2 = f(x + h/2*k1);
    k3 = f(x + h/2*k2);
    k4 = f(x + h*k3);
    x_next = x + h/6*(k1 + 2*k2 + 2*k3 + k4);
end

